#!/usr/bin/env python3
from __future__ import division
import sys
import os
import time
import random
import re
import subprocess

from Bio.Blast import NCBIXML
from Bio import SeqIO

import collections


class Blaster():
    def __init__(self, inputfile, databases, db_path, out_path='', min_cov=0.6,
                 threshold=0.9, blast='blastn', cut_off=True,
                 max_target_seqs=50000, reuse_results=False,
                 allowed_overlap=0):

        min_cov = 100 * float(min_cov)
        threshold = 100 * float(threshold)

        # For alignment data storage
        self.gene_align_query = dict()  # Sequence alignment lines
        self.gene_align_homo = dict()  # Sequence alignment homolog string
        self.gene_align_sbjct = dict()  # Sequence alignment allele string
        self.results = dict()  # Results

        # TODO: Add excluded results to this dictionay
        self.results["excluded"] = dict()

        for db in databases:
            # Adding the path to the database and output
            db_file = "%s/%s.fsa" % (db_path, db)
            tmp_out_path = "%s/tmp" % (out_path)
            out_file = "%s/out_%s.xml" % (tmp_out_path, db)
            os.makedirs(tmp_out_path, exist_ok=True)
            os.chmod(tmp_out_path, 0o775)

            # Running blast
            if (os.path.isfile(out_file)
                    and os.access(out_file, os.R_OK)
                    and reuse_results):
                print("Found " + out_file + " skipping DB.")
                out, err = (b'', b'')
            else:
                cmd = ("%s -subject %s -query %s -out %s -outfmt '5'"
                       " -perc_identity  %s -max_target_seqs %s"
                       " -dust 'no'" % (blast, db_file, inputfile,
                                        out_file, threshold, max_target_seqs))

                process = subprocess.Popen(cmd, shell=True,
                                           stdout=subprocess.PIPE,
                                           stderr=subprocess.PIPE)
                out, err = process.communicate()

            # Get results file
            try:
                 # Test if output file exist
                 result_handle = open(out_file, "r")
            except IOError:
                 sys.exit(("Error: BLAST did not run as expected. "
                           "The expected output file, {}, did not exist.\n"
                           "BLAST finished with the following response:"
                           "\n{}\n{}").format(os.path.abspath(out_file),
                                              out.decode("utf-8"),
                                              err.decode("utf-8")))

            # Test if blast output is empty
            if os.stat(out_file).st_size == 0:
                 sys.exit(("Error: BLAST did not run as expected. "
                           "The output file {} was empty.\n"
                           "BLAST finished with the following response:"
                           "\n{}\n{}").format(os.path.abspath(out_file),
                                              out.decode("utf-8"),
                                              err.decode("utf-8")))

            # Get blast output
            blast_records = NCBIXML.parse(result_handle)

            # Declaring variables for saving the results

            # results for each gene
            gene_results = dict()

            # For finding the best hits
            best_hsp = dict()

            # Keeping track of gene split
            gene_split = collections.defaultdict(dict)

            # Making the dicts for sequence outputs
            self.gene_align_query[db] = dict()
            self.gene_align_homo[db] = dict()
            self.gene_align_sbjct[db] = dict()

            # Parsing over the hits and only keeping the best
            for blast_record in blast_records:
                # OLD CODE TO BE REMOVED
                #  query = blast_record.query
                #  blast_record.alignments.sort(key=lambda align: (
                #     -max((len(hsp.query)
                #            * (int(hsp.identities) / float(len(hsp.query)))
                #            for hsp in align.hsps)))
                #  )

                # Sort BLAST alignments by the hsp in each alignment with the
                # highest number of identical nucleotides (hsp.identities)
                blast_record.alignments.sort(
                    key=lambda align: (max((int(hsp.identities)
                                           for hsp in align.hsps))),
                    reverse=True)
                query = blast_record.query

                for alignment in blast_record.alignments:
                    # Setting the e-value as 1 and bit as 0 to get the best
                    # HSP fragment
                    best_e_value = 1
                    best_bit = 0
                    start_hsp = 0
                    end_hsp = 0

                    for hsp in alignment.hsps:
                        if hsp.expect < best_e_value or hsp.bits > best_bit:
                            # best_e_value = hsp.expect
                            # best_bit = hsp.bits
                            tmp = alignment.title.split(" ")
                            sbjct_header = tmp[1]
                            # DEBUG
                            print("Found: {}".format(sbjct_header))
                            bit = hsp.bits
                            sbjct_length = alignment.length
                            sbjct_start = hsp.sbjct_start
                            sbjct_end = hsp.sbjct_end
                            gaps = hsp.gaps
                            query_string = str(hsp.query)
                            homo_string = str(hsp.match)
                            sbjct_string = str(hsp.sbjct)
                            contig_name = query.replace(">", "")
                            query_start = hsp.query_start
                            query_end = hsp.query_end
                            HSP_length = len(query_string)
                            perc_ident = (int(hsp.identities)
                                          / float(HSP_length) * 100)
                            strand = 0
                            coverage = ((int(HSP_length) - int(gaps))
                                        / float(sbjct_length))
                            perc_coverage = (((int(HSP_length) - int(gaps))
                                             / float(sbjct_length)) * 100)

                            # cal_score is later used to select the best hit
                            cal_score = perc_ident * coverage

                            hit_id = "%s:%s..%s:%s:%f" % (
                                contig_name, query_start, query_end,
                                sbjct_header, cal_score)
                            # If the hit is on the other strand
                            if sbjct_start > sbjct_end:
                                tmp = sbjct_start
                                sbjct_start = sbjct_end
                                sbjct_end = tmp

                                query_string = self.reversecomplement(
                                    query_string)
                                homo_string = homo_string[::-1]
                                sbjct_string = self.reversecomplement(
                                    sbjct_string)
                                strand = 1

                            # Save hit
                            if((cut_off and perc_coverage > 20)
                               or cut_off is False):

                                best_hsp = {'evalue': hsp.expect,
                                            'sbjct_header': sbjct_header,
                                            'bit': bit,
                                            'perc_ident': perc_ident,
                                            'sbjct_length': sbjct_length,
                                            'sbjct_start': sbjct_start,
                                            'sbjct_end': sbjct_end,
                                            'gaps': gaps,
                                            'query_string': query_string,
                                            'homo_string': homo_string,
                                            'sbjct_string': sbjct_string,
                                            'contig_name': contig_name,
                                            'query_start': query_start,
                                            'query_end': query_end,
                                            'HSP_length': HSP_length,
                                            'coverage': coverage,
                                            'cal_score': cal_score,
                                            'hit_id': hit_id,
                                            'strand': strand,
                                            'perc_coverage': perc_coverage
                                            }



                        # Saving the result if any
                        if best_hsp:
                            save = 1

                            # If there are other gene alignments they are compared
                            if gene_results:
                                tmp_gene_split = gene_split
                                tmp_results = gene_results
                                # Compare the hit results
                                save, gene_split, gene_results = (
                                    self.compare_results(save, best_hsp,
                                                         tmp_results,
                                                         tmp_gene_split,
                                                         allowed_overlap)
                                )

                            # If the hit is not overlapping with other hit
                            # seqeunces it is kept
                            if save == 1:
                                # DEBUG
                                print("Saving: {}".format(hit_id))
                                gene_results[hit_id] = best_hsp

            result_handle.close()

            # If the hit does not cover the entire database reference the
            # missing seqence data are extracted
            keys = list(gene_results.keys())

            for hit_id in keys:
                hit = gene_results[hit_id]

                # Calculate possible split gene coverage
                perc_coverage = hit['perc_coverage']

                if(hit['sbjct_header'] in gene_split
                    and len(gene_split[hit['sbjct_header']]) > 1):

                    # Calculate new length
                    new_length = self.calculate_new_length(gene_split, gene_results,
                                                                        hit)
                    hit['split_length'] = new_length

                    # Calculate new coverage
                    perc_coverage = new_length / float(hit['sbjct_length']) * 100

                # If the hit is above the minimum length threshold it is kept
                if perc_coverage >= min_cov:
                    if hit['coverage'] == 1:
                        self.gene_align_query[db][hit_id] = hit['query_string']
                        self.gene_align_homo[db][hit_id] = hit['homo_string']
                        self.gene_align_sbjct[db][hit_id] = hit['sbjct_string']
                    elif hit['coverage'] != 1:
                        # Getting the whole database sequence
                        for seq_record in SeqIO.parse(db_file, "fasta"):
                            if seq_record.description.replace(" ", "") == hit['sbjct_header'].replace(" ", ""):
                                start_seq = str(seq_record.seq)[:int(hit["sbjct_start"])-1]
                                end_seq = str(seq_record.seq)[int(hit["sbjct_end"]):]
                                self.gene_align_sbjct[db][hit_id] = start_seq + hit['sbjct_string'] + end_seq
                                #self.gene_align_sbjct[db][hit_id] = str(seq_record.seq)
                                break

                        # Getting the whole contig to extract extra query seqeunce
                        contig = ''
                        for seq_record in SeqIO.parse(inputfile, "fasta"):
                            if seq_record.description.replace(" ", "") == hit['contig_name'].replace(" ", ""):
                                contig = str(seq_record.seq)
                                break

                        # Extract extra sequence from query
                        query_seq, homo_seq = self.get_query_align(hit, contig)

                        # Saving the new alignment sequences
                        self.gene_align_query[db][hit_id] = query_seq
                        self.gene_align_homo[db][hit_id] = homo_seq

                else:
                    del gene_results[hit_id]
                    if hit['sbjct_header'] in gene_split:
                        del gene_split[hit['sbjct_header']]

            # Save the database result
            if gene_results:
                self.results[db] = gene_results
            else:
                self.results[db] = "No hit found"

    @staticmethod
    def reversecomplement(seq):
        # Make reverse complement strand
        trans = str.maketrans("ATGC", "TACG")
        return seq.translate(trans)[::-1]

    @staticmethod
    def compare_results(save, best_hsp, tmp_results, tmp_gene_split,
                        allowed_overlap):
        """
            Function for comparing hits and saving only the best hit
        """

        # Get data for comparison
        hit_id = best_hsp['hit_id']
        new_start_query = best_hsp['query_start']
        new_end_query = best_hsp['query_end']
        new_start_sbjct = int(best_hsp['sbjct_start'])
        new_end_sbjct = int(best_hsp['sbjct_end'])
        new_score = best_hsp['cal_score']
        new_db_hit = best_hsp['sbjct_header']
        new_contig = best_hsp['contig_name']
        new_HSP = best_hsp['HSP_length']

        # See if the best HSP fragment overlap with another allignment
        # and keep the allignment with the highest score - if the new
        # fragment is not providing new sequence
        keys = list(tmp_results.keys())
        for hit in keys:
            hit_data = tmp_results[hit]
            old_start_query = hit_data['query_start']
            old_end_query = hit_data['query_end']
            old_start_sbjct = int(hit_data['sbjct_start'])
            old_end_sbjct = int(hit_data['sbjct_end'])
            old_score = hit_data['cal_score']
            old_db_hit = hit_data['sbjct_header']
            old_contig = hit_data['contig_name']
            old_HSP = hit_data['HSP_length']

            remove_old = 0
            # If they align to the same gene in the database they are
            # compared
            if new_db_hit == old_db_hit:

                #If the hit comes from different contig
                if old_contig != new_contig:

                    # Save a split if the new hit still creats one
                    if(new_db_hit in tmp_gene_split
                       and hit_id not in tmp_gene_split[new_db_hit]):

                        tmp_gene_split[new_db_hit][hit_id] = 1

                # If the hit provides additional sequence it is kept and
                # the new coverage is saved otherwise the one with the
                # highest score is kept
                elif(new_start_sbjct < old_start_sbjct
                   or new_end_sbjct > old_end_sbjct):

                    # Save the hits as split
                    tmp_gene_split[old_db_hit][hit_id] = 1
                    if hit not in tmp_gene_split[old_db_hit]:
                        tmp_gene_split[old_db_hit][hit] = 1


#                else:
#                    if new_score > old_score:
                        # Set to remove old hit
#                        remove_old = 1

                        # Save a split if the new hit still creats one
#                        if(new_db_hit in tmp_gene_split
#                           and hit_id not in tmp_gene_split[new_db_hit]):

#                            tmp_gene_split[new_db_hit][hit_id] = 1
#                    else:
#                        save = 0

                        # If the old and new hit is not identical the
                        # possible saved gene split for the new hit is
                        # removed
#                        if hit_id != hit:
#                            if(new_db_hit in tmp_gene_split
#                               and hit_id in tmp_gene_split[new_db_hit]):

#                                del tmp_gene_split[new_db_hit][hit_id]
#                        break

            # If the hits comes form the same part of the contig
            # sequnce but match different genes only the best hit is
            # kept
            if new_contig == old_contig:
                print("Same contig: {} == {}".format(new_contig, old_contig))
                print("\t{} vs {}".format(new_db_hit, old_db_hit))
                # Check if saved hits overlaps with current hit
                hit_union_length = (max(old_end_query, new_end_query)
                                    - min(old_start_query, new_start_query))
                hit_lengths_sum = ((old_end_query - old_start_query)
                                   + (new_end_query - new_start_query))
                overlap_len = (hit_lengths_sum - hit_union_length)

                if overlap_len < allowed_overlap:
                    print("\tignore overlap ({}): {}".format(overlap_len, new_db_hit))
                    continue
                print("\toverlap found ({}): {}".format(overlap_len, new_db_hit))

                # If the two hits cover the exact same place on the
                # contig only the percentage of identity is compared
                if(old_start_query == new_start_query
                   and old_end_query == new_end_query):

                    if best_hsp['perc_ident'] > hit_data['perc_ident']:

                        # Set to remove old hit
                        remove_old = 1

                        # Save a split if the new hit still creats one
                        if(new_db_hit in tmp_gene_split
                           and hit_id not in tmp_gene_split[new_db_hit]):

                            tmp_gene_split[new_db_hit][hit_id] = 1

                    elif best_hsp['perc_ident'] == hit_data['perc_ident']:
                        # Save both

                        # Save a split if the new hit still creats one
                        if(new_db_hit in tmp_gene_split
                           and hit_id not in tmp_gene_split[new_db_hit]):

                            tmp_gene_split[new_db_hit][hit_id] = 1
                    else:
                        save = 0

                        # Remove new gene from gene split if present
                        if(new_db_hit in tmp_gene_split
                           and hit_id in tmp_gene_split[new_db_hit]):

                            del tmp_gene_split[new_db_hit][hit_id]
                        break

                # If new hit overlaps with the saved hit
                elif(hit_union_length <= hit_lengths_sum):
                    print("\t{} <= {}".format(hit_union_length, hit_lengths_sum))
                    print("\t\tScores: {} cmp {}".format(new_score, old_score))

                    if new_score > old_score:
                        # Set to remove old gene
                        remove_old = 1
                        # Save a split if the new hit still creats one
                        if(new_db_hit in tmp_gene_split
                           and hit_id not in tmp_gene_split[new_db_hit]):

                            tmp_gene_split[new_db_hit][hit_id] = 1

                    elif new_score == old_score:
                        # If both genes are of same coverage
                        # and identity is the same implied by new_score == old_score
                        # hit is chosen based on length
                        if((int(best_hsp['perc_coverage']) ==
                             int(hit_data['perc_coverage']))
                             and new_HSP > old_HSP):

                            # Set to remove old gene
                            remove_old = 1
                        elif((int(best_hsp['perc_coverage']) ==
                                int(hit_data['perc_coverage']))
                                and old_HSP > new_HSP):

                            # Remove current hit
                            save = 0
                        elif((int(best_hsp['perc_coverage']) ==
                                int(hit_data['perc_coverage']))
                                and old_HSP==new_HSP):
                            # Both hits has same coverage, and same identity
                            # and same length, how to choose only one hit?
                            pass

                        # TODO
                        # If new_score == old_score but identity and coverages are not the same.
                        # which gene should be chosen?? Now they are both keept.

                        # Save a split if the new hit creats one - both
                        # hits are saved
                        if(new_db_hit in tmp_gene_split
                            and hit_id not in tmp_gene_split[new_db_hit]):

                            tmp_gene_split[new_db_hit][hit_id] = 1
                    else:
                        # Remove new gene from gene split if present
                        if(new_db_hit in tmp_gene_split
                            and hit_id in tmp_gene_split[new_db_hit]):

                            del tmp_gene_split[new_db_hit][hit_id]

                        save = 0
                        break

            # Remove old hit if new hit is better
            if remove_old == 1:
                del tmp_results[hit]

                # Remove gene from gene split if present
                if(old_db_hit in tmp_gene_split
                    and hit in tmp_gene_split[old_db_hit]):

                    del tmp_gene_split[old_db_hit][hit]

        return save, tmp_gene_split, tmp_results

    @staticmethod
    def calculate_new_length(gene_split, gene_results, hit):
        """
            Function for calcualting new length if the gene is split on
            several contigs
        """
        # Looping over splitted hits and calculate new length
        first = 1
        for split in gene_split[hit['sbjct_header']]:

            new_start = int(gene_results[split]['sbjct_start'])
            new_end = int(gene_results[split]['sbjct_end'])

            # Get the first HSP
            if first == 1:
                new_length = int(gene_results[split]['HSP_length'])
                old_start = new_start
                old_end = new_end
                first = 0
                continue
            if new_start < old_start:
                new_length = new_length + (old_start - new_start)
                old_start = new_start

            if new_end > old_end:
                new_length = new_length + (new_end - old_end)
                old_end = new_end

        return(new_length)

    @staticmethod
    def get_query_align(hit, contig):
        """
            Function for extracting extra seqeunce data to the query
            alignment if the full reference length are not covered
        """

        # Getting data needed to extract sequences
        query_seq = hit['query_string']
        homo_seq = hit['homo_string']
        sbjct_start = int(hit['sbjct_start'])
        sbjct_end = int(hit['sbjct_end'])
        query_start = int(hit['query_start'])
        query_end = int(hit['query_end'])
        length = int(hit['sbjct_length'])

        # If the alignment doesn't start at the first position data is
        # added to the begnning
        if sbjct_start != 1:
            missing = sbjct_start - 1

            if(query_start >= missing and hit['strand'] != 1
                or hit['strand'] == 1 and missing <= (len(contig) - query_end)):

                # Getting the query sequence.
                # If the the hit is on the other strand the characters
                # are reversed.
                if hit['strand'] == 1:
                    start_pos = query_end
                    end_pos = query_end + missing
                    chars = contig[start_pos:end_pos]
                    chars = Blaster.reversecomplement(chars)
                else:
                    start_pos = query_start - missing - 1
                    end_pos = query_start - 1
                    chars = contig[start_pos:end_pos]

                query_seq = chars + str(query_seq)
            else:
                # Getting the query sequence.
                # If the the hit is on the other strand the characters
                # are reversed.
                if hit['strand'] == 1:
                    if query_end == len(contig):
                        query_seq = "-" * missing + str(query_seq)
                    else:
                        start_pos = query_end
                        chars = contig[start_pos:]
                        chars = Blaster.reversecomplement(chars)

                        query_seq = ("-" * (missing - len(chars))
                                         + chars + str(query_seq))
                elif query_start < 3:
                    query_seq = "-" * missing + str(query_seq)
                else:
                    end_pos = query_start - 2
                    chars = contig[0:end_pos]

                    query_seq = ("-" * (missing - len(chars))
                                     + chars + str(query_seq))

            # Adding to the homo sequence
            spaces = " " * missing
            homo_seq = str(spaces) + str(homo_seq)

        # If the alignment dosen't end and the last position data is
        # added to the end
        if sbjct_end < length:
            missing = length - sbjct_end

            if(missing <= (len(contig) - query_end) and hit['strand'] != 1
                or hit['strand'] == 1 and query_start >= missing):

                # Getting the query sequence.
                # If the the hit is on the other strand the characters
                # are reversed.
                if hit['strand'] == 1:
                    start_pos = query_start - missing - 1
                    end_pos = query_start - 1
                    chars = contig[start_pos:end_pos]
                    chars = Blaster.reversecomplement(chars)
                else:
                    start_pos = query_end
                    end_pos = query_end + missing
                    chars = contig[start_pos:end_pos]

                query_seq = query_seq + chars
            else:
                # If the hit is on the other strand the characters are reversed
                if hit['strand'] == 1:
                    if query_start < 3:
                        query_seq = query_seq + "-" * missing
                    else:
                        end_pos = query_start - 2
                        chars = contig[0:end_pos]
                        chars = Blaster.reversecomplement(chars)

                        query_seq = (query_seq
                                         + chars + "-" * (missing - len(chars)))
                elif query_end == len(contig):
                    query_seq = query_seq + "-" * missing
                else:
                    start_pos = query_end
                    chars = contig[start_pos:]

                    query_seq = query_seq + chars + "-" * (missing - len(chars))

            # Adding to the homo sequence
            spaces = " " * int(missing)
            homo_seq = str(homo_seq) + str(spaces)

        return query_seq, homo_seq
